Spontaneous symmetry breaking in a generalized orbital compass model 
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We introduce a generalized two-dimensional orbital compass model, which interpolates contin- 
uously from the classical Ising model to the orbital compass model with frustrated quantum in- 
teractions, and investigate it using the multiscale entanglement renormalization ansatz (MERA). 
The results demonstrate that increasing frustration of exchange interactions triggers a second or- 
der quantum phase transition to a degenerate symmetry broken state which minimizes one of the 
interactions in the orbital compass model. Using boson expansion within the spin-wave theory we 
unravel the physical mechanism of the symmetry breaking transition as promoted by weak quan- 
tum fluctuations and explain why this transition occurs only surprisingly close to the maximally 
frustrated interactions of the orbital compass model. The spin waves remain gapful at the critical 
point, and both the boson expansion and MERA do not find any algebraically decaying spin- spin 
correlations in the critical ground state. 
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I. INTRODUCTION 

The orbital compass model (OCM) is physically moti- 
vated by the orbital interactions which arise for strongly 
correlated electrons in transition metal oxides with partly 
filled degenerate 3d orbitals and lead to rich and still 
poorly understood quantum models. In these systems 
the orbital degrees of freedom play a crucial role in deter- 
mining collective states such as coexisting magnetic and 
orbital order, as for instance in the colossal magnetoresis- 
tance manganites 1 or in the vanadate perovskites^ The 
orbital interactions are typically intrinsically frustrated 
and may strongly enhance quantum fluctuations, leading 
to disordered states^ While realistic orbital interactions 
are somewhat complex, a paradigm of intrinsic frustra- 
tion is best realized in the OCM, 4 8 with the pseudospin 
couplings intertwined with the orientation of interacting 
bonds. Its two-dimensional (2D) version on a honeycomb 
lattice^ realized in layered iron oxides, 10 is equivalent to 
the Kitaev model. 11 

Although conceptually quite simple, the OCM has an 
interdisciplinary character as it plays an important role 
in a variety of contexts beyond the correlated transi- 
tion metal oxides, such as: (i) the implementation of 
protected qubits for quantum computations in Joseph- 
son lattice arrays j£ (ii) topological quantum order ji 2 - or 
(iii) polar molecules in optical lattices and systems of 
trapped ions* 1 ^ Numerical studies 7 suggested that when 
anisotropic interactions are varied through the isotropic 
point of the 2D OCM, the ground state is not an or- 
bital liquid type but instead a first order quantum phase 
transition (QPT) occurs between two different types of 
Ising-type order dictated by one or the other interac- 
tion. Recently the existence of this transition, similar 
to the one which occurs in the exact solution of the 
one-dimensional OCM, 14 was confirmed using projected 



entangled-pair state algorithm. 15 This implies that the 
symmetry is spontaneously broken at the compass point, 
and the spin order follows one of the two equivalent frus- 
trated interactions. 

Knowing that the ground states of the 2D Ising model 
and the 2D OCM are quite different, we introduce a gen- 
eralized OCM which interpolates between these two lim- 
iting cases. Using this model we will investigate: (i) 
the physical consequences of gradually increasing frus- 
tration in a 2D system, (ii) where a QPT occurs from 
the Ising ground state to the degenerate ground state 
of the OCM, and, finally, (iii) the order and the physi- 
cal mechanism of this QPT. As increasing frustration of 
the orbital interactions introduces entangled states, the 
present problem provides a unique opportunity to use the 
recently developed multiscale entanglement renormaliza- 
tion ansatz 16,17 (MERA) in order to find reliable answers 
to the above questions. As we show below, the QPT in 
the generalized OCM occurs only surprisingly close to 
the maximally frustrated interactions in the OCM. We 
also explain the physical origin of this behavior using an 
analytic approach based on the spin- wave theory. 

Quantum many-body systems exhibit several interest- 
ing collective phenomena. Recent progress in develop- 
ing efficient numerical methods to study quantum sys- 
tems on a lattice is remarkable and allowed to investi- 
gate complex many-body phenomena, including QPTs^ 
An important step here was the discovery of density ma- 
trix renormalization group, 19 a very powerful numerical 
method that can be applied to one-dimensional strongly 
correlated fermionic and bosonic systems. 20 This idea 
played a fundamental role in developing entanglement 
renormalizatio n 16 ! 17 to study quantum spin systems on 
a 2D lattice. Crucial in this approach is the removal 
of short-range entanglement by unitary transformations 
called disent anglers. It generates a real-space renormal- 
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ization group transformation implemented in the MERA, 
which was recently successfully employed to investi- 
gate several quantum spin models^"— and interacting 
fermion systems i 26 i 27 So far, the very promising MERA 
has been applied inter alia to the 2D quantum Ising 
model s 21 ! 22 and to the Heisenberg model on a kagome 
lattice, 28 but other possible applications and the optimal 
geometries for performing sequentially disentanglement 
and isometry transformation were also discussed^ 

The paper is organized as follows. In Sec. II we in- 
troduce the generalized OCM and state the problem of 
the existence and nature of the QPT. Next we present 
the MERA algorithm in Sec. IHII used to investigate frus- 
trated interactions in the OCM. Numerical results ob- 
tained using the MERA are presented in Sec. [IVJ In or- 
der to explain the physical mechanism of the QPT found 
in the OCM we performed the boson expansion within 
the spin- wave theory, as described in Sec. [Vj The details 
of this expansion are presented in the Appendix. The 
paper is summarized in Sec. IVI[ where the main conclu- 
sions of the present work are also given. 



II. GENERALIZED COMPASS MODEL 

In this paper, we investigate the nature and position 
of the QPT when the OCM point is approached in a 
different way from that studied before, namely when 
frustration of interactions along two nonequivalent di- 
rections gradually increases. Therefore, we introduce a 
2D generalized OCM with ferro-like interactions 29 on a 
square lattice in ab plane (we assume the exchange con- 
stant J = 1), 

no) = -£ {<-(0K+u (0) + (*)} • (!) 

ijGab 

The interactions occur between nearest neighbors and 
are balanced along both lattice directions a and b. Here 
{ij} labels lattice sites, with i (j) increasing along a (b) 
axes, and {crfj(6), cr^.(0)} are linear combinations of Pauli 
matrices describing interactions for S = 1/2 spins: 

atj (0) = cos(0/2) af s + sin(0/2) afj , (2) 
a% (0) = cos(0/2) of j - sin(0/2) <r?- . (3) 

The interactions in Eq. (pQ) include the classical Ising 
model at 6 = 0° for afj operators and become gradu- 
ally more frustrated with increasing angle G (0°,90°] 
— they interpolate between the Ising model (at = 0°) 
and the isotropic OCM (at = 90°), see Fig. [Q The 
latter case is equivalent to the 2D OCM with standard 
interactions af-af-^ and afjaf +1 j along the a and b 
directions 4 8 by a straightforward unitary transforma- 
tion. The model (pQ) includes also as a special case the 
2D orbital model for e g electrons at = 60°, 29 describ- 
ing, for instance, the orbital part of the superexchange 
interactions in the ferromagnetic planes of LaMn03^S 



(a) = 0° (b) = 60° (c) = 90° 




FIG. 1: (Color online) Artist's view of the evolution of orbital 
interactions in the generalized OCM Eq. (pQ) with increasing 
angle 0. Heavy (blue) lines indicate favored spin direction 
induced by interactions along two nonequivalent lattice axes 
a and b. Different panels show: (a) the Ising model at = 0°, 
(b) the 2D e g orbital model at = 60°, and (c) the OCM at 
— 90°. Spin order follows the interactions in the Ising limit, 
while it follows one of the equivalent interactions, a a or a b , in 
the OCM. This results in the symmetry breaking QPT which 
occurs between (b) and (c), as we show in Sees. llVfVl 



Since the isotropic model has the same interaction 
strength for the bonds along both a and b axis, it is sym- 
metric under transformation a o 6, and the issue of the 
QPT between different ground states of the anisotropic 
compass model 15 does not arise. On one hand, this sym- 
metry is obeyed by the classical Ising ground state, while 
on the other hand, in the ground state of the OCM this 
symmetry is spontaneously broken (and the ground state 
is degenerate). Therefore, an intriguing question con- 
cerning the ground state of the model (pp) is whether it 
has the same high symmetry as the Ising model in a broad 
range of or the symmetry is soon spontaneously broken 
when increases, i.e., there are degenerate ground states 
with lower symmetries, also for the e g orbital model, see 
Fig. mb). This question has been addressed by inves- 
tigating the energy contributions along two equivalent 
lattice directions a and b by applying the MERA. 



III. MERA ALGORITHM 

A. Calculation method 

In order to obtain the ground state, we use a transla- 
tionally invariant MERA on infinite lattice. 23 The MERA 
is a tensor network with infinite number of layers of disen- 
t anglers and isometries. By translational invariance, all 
isometries (disent anglers) in a given layer are the same. 
Since every layer represents a coarse-graining renormal- 
ization group transformation, shown in Fig. [2] for the 
9-to-l geometry and described in more detail in Ref. [23| 
(see their Fig. 7), we assume that after a finite number of 
such transformations a fixed point of the renormalization 
group is reached (either trivial or non-trivial) and from 
that time on the following transformations are the same. 
In other words, at the bottom of the tensor network there 
is a finite number of non-universal layers whose tensors 
are different in general, but above certain level all layers 
are the same. The bottom layers describe non-universal 
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FIG. 2: (Color online) 9-to-l geometry of MERA applied to 
the OCM: dark (red) boxes represent the action of disent an- 
glers U, Uh, U v and gray (green) ones — isometries W; arrows 
indicate subsequent transformations used; the labels of spins 
1-9 in a single block are addressed in the text. This is a 
coarse-graining renormalization group transformation where 
each 3x3 plaquette in the top-left panel is replaced by a 
coarse grained spin in the bottom-right panel (nine spins are 
replaced by one coarse-grained spin). To minimize the num- 
ber of states x of the coarse-grained spin, the microscopic 
spins are disentangled prior to decimation. We also used a 
5-to-l geometry, see Fig. 1 of Ref. I22I . 



short range correlations, and the universal layers above 
this level describe universal properties of the fixed point. 
The number N of the non-universal bottom layers is one 
of the parameters of the infinite-lattice MERA. We have 
verified that it is enough to keep up to three non-universal 
layers, depending on how close the critical point is. 

Starting with randomly chosen tensors, the structure 
is optimized layer by layer, from the top to bottom and 
back. In given layer r, we calculate an environment of 
each tensor type by means of renormalized Hamiltonians 
h T and density matrices p T computed from other lay- 
ers. The environments are aimed at updating tensors to 
minimize total energy. In the universal layer, this updat- 
ing technique is slightly different: and are fixed 
points of the renormalization procedure defined by ten- 
sors in this layer. The above steps are iterated until the 
convergence of energy is achieved. For given we obtain 
the ground state for different values of bond dimension 
X- It turns out that in most cases it is sufficient to work 
with x = 3, which is the same in each layer. However, 
it is necessary to increase x to 4 in the neighborhood 
of the critical point. The number of operations and the 
required memory scale as x 16 and x 12 respectively. The 
inset in Fig. H(a) shows the convergence of the energy 
of the ground state with an increasing bond dimension. 
Here we also present a comparison of results obtained 
with the alternative 5-to-l geometry. 

The algorithm is implemented in C++ and optimized 
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FIG. 3: (Color online) (a) Renormalizer superoperator which 
consists of isometries only. The connections show how the 
isometries are contracted; compare Eq. (0. (b) Method of 
calculating correlations (o x o y ) = (00O0) between two sites 
separated by a distance 3 n+iV . The scheme presents a graph- 
ical explanation of Eq. ©. (c) Deriving p& from />oo when 
sites are separated vertically (top) and horizontally (bottom) . 



in order to work on multi-processor computers. On an 
eight-core 2.3 GHz processor, it takes about half an hour 
to update the whole tensor network which consists of four 
layers of tensors with x = 4. Near the critical point, i.e., 
at ~ # c , the convergence requires several thousands 
of iterations whereas it is significantly faster far from C . 
When is scanned from 0° to 90° (or back), it is more ef- 
ficient to use the previous ground state as an initial state 
for the next discrete value of instead of starting from 
a random initial state for each value of 6. We have care- 
fully verified convergence to the ground state by scanning 
back and forth and comparing the results with those 
obtained from random initial states for selected values of 



B. Correlations 

In order to calculate correlations, we take advantage of 
the special structure of the renormalization group trans- 
formation in Fig. [21 A site of the lattice that lies in 
the center of a 3 x 3 decimation block (number 9 in Fig. 
[2]) undergoes renormalization in a particularly easy man- 
ner. Since no disent angler is applied to this central site, a 
one-site operator o r _i at this site is mapped by the r-th 
renormalization group transformation to a coarse-grained 
one-site operator 



R T o T - 



(4) 



Here R T is a renormalizer superoperator built out of con- 
tracted isometries only, as shown in Fig. [3fa): 



( R r)kl 



m ...ngl 



(5) 



The meaning of the transformations W T and W} is given 
in Figs. [2] and [3](a). Thus, if we have N non- universal 
layers at the bottom of the geometry of MERA, then 
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renormalized one-site operators at the central sites just 
below the universal layer are given by [see Fig. [3fb)]: 



on 



R N R 



N-l 



(6) 



where oo = o denotes a physical, microscopic one-site 
operator at one of the central sites at the very bottom of 
the MERA tensor network. 

To extract information on the correlations, it is con- 
venient to write eigen-decomposition of the renormalizer 
Roo in the universal layer: 



Ro 



(7) 



It is straightforward to verify the basic property of the 
spectrum of R^: \X a \ < 1. The ortonormality of the 
vectors W l in Eq. (j5]) implies that the identity operator 
( v i)ij — $ij is an eigenvector with eigenvalue Ai = 1. 
In our numerical calculations this is the only eigenvalue 
with modulus 1. 

After the operator on is decomposed as on = 
^2a°N v ai a repeated action of the renormalizer in 
the universal layers can be written as 



R^qOn 



°N v a 



(8) 



A correlator between two central sites x and y separated 
by a distance |x — y| = 3 n+Ar in the horizontal (vertical) 
direction is thus given by: 

<o x o y > = Tr (R^o N R^o N )} (9) 
= J2°>N^X n a \ n (10) 

a, (3 



where 



3 n and 



at,/3 



Tr 



a P 
On°N C aP 
r - log 3 (A a A /3 ) ' 



(11) 



(12) 



Here poo is a two-site reduced density matrix in a uni- 
versal layer derived from as depicted in Fig. [3fc). 

Correlations corresponding to the leading eigenvalue 
Ai = 1 do not decay with the distance between x and y. 
They describe long range order in the operator o and can 
be used to extract its expectation value (o): 

(o) 2 = lim (o x o y ) = o^o^cn = (oj^) 2 , (13) 

|x-y|^oo 

where we use the property: lim n ^oo = that holds for 
a > 1 and the fact that en = 1 which is a consequence 
of vi being an identity. Thus only a one-site operator 
with a non-zero coefficient o]^ has non-zero expectation 
value. A trivial example is the identity o = I. Indeed, we 
obtain oat = I in Eq. (J6j), which is equivalent to o]^ — 1, 
and Eq. (fT3|) yields (I) 2 = 1 as expected. 



IV. NUMERICAL RESULTS 
A. Symmetry breaking transition 

Information about the ground state of the OCM Eq. 
(PQ) is contained in average energy per bond E(0) and 
energy anisotropy AE(0): 



E ( ) = -o <^K+ijW + <W<i + iW)>(14) 



AE(0) 



«(0K + i,#)} - (4( 6 Hj + i(°)) -( 15 ) 



In the classical limit of Ising interactions E(0°) = —1 and 
A£?(0°) = 0. Due to increasing frustration, the energy 
E{0) gradually increases for increasing angle 9 in Eq. (pp) 
and reaches a maximum of £?(90°) ~ —0.57 in the OCM, 
see Fig. H^a). This increase is smooth and does not 
indicate the existence of a QPT. 

However, by investigating the anisotropy AE(6) Eq. 
(fT5j) between a and b bonds, we identified an angle 6 C at 
which AE(0) starts to grow. Although a gradual evolu- 
tion of the ground state staring from = 0° might be also 
expected, the Ising-type state is first surprisingly robust 
in a broad range of angles G [O°,0 C ], and the energy 
associated with bonds along the a and b axes remains 
the same, i.e., AE(0) = 0. Next, the symmetry between 
the a and b directions is spontaneously broken above # c , 
where a finite value of AE{0) is found, and then AE{9) 
grows rapidly with further increasing angle i.e., large 
spin correlations develop along only one of the two equiv- 
alent directions a and b. This QPT was detected by the 
MERA at 9 C ~ 84.8°, see Fig. g^b). 



B. Magnetization in the ground state 

To understand better the QPT at C let us consider 
the expectation value of the spontaneous magnetization 
M = {M x , M y , M z } derived from the long range order 
in the correlation function: 



lim 

|x-yK 



(<j*(T l y ) = M k M l 



(16) 



where k(l) = x,y,z. For the interactions in Eq. (1) one 
finds M y = for any 0. 

We found that the ground state obtained using the 
MERA for < C is characterized by M z = and Ising- 
type long range order of M x which gradually decreases 
but remains rather large, \M X \ > 0.93, in this parameter 
range. The symmetry between the directions a and b is 
broken above C by appearance of a nonzero component 
M z . 

The value of the total magnetization 



M = |M| = \J (M x ) 2 + (M z ) 2 , 



(17) 



obtained from the MERA decreases continuously from 
M (0°) = 1 in the Ising model to M(90°) ~ 0.92 in the 
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FIG. 4: (Color online) Ground state obtained for the gener- 
alized OCM Eq. (1) using the MERA: (a) average energy E 
per bond given by Eq. (HJ), (b) energy anisotropy AE given 
by Eq. (|15p , (c) spontaneous magnetization M given by Eq. 
(|17p . and (d) magnetization orientation <f> given by Eq. (|18[). 
Embedded L x L clusters coupled to the neighboring spins 
by mean-field terms (L x L MF) exhibit qualitatively simi- 
lar behavior. Inset: Convergence of the ground state energy 
obtained by two geometries of MERA with increasing bond 
dimension x- Black: 9-to-l geometry presented in Fig. [2j 
blue: 5-to-l geometry introduced in Fig. 1 of Ref. [22|. The 
9-to-l geometry results prove to converge faster for close to 



OCM, see Fig. HJ^c). Thus the reduction in the order 
parameter M by quantum fluctuations arising from the 
admixture of the z-th component, is here rather small, 
and reproduces qualitative results obtained for the e g or- 
bital model within the linear orbital wave theory^!. Fur- 
thermore, by a closer inspection of M{0) we have found 
that the derivative (8M(0)/dO) does not exist at = C . 

As expected from the behavior of AE, the obtained 
symmetry breaking shown in Fig. [4] implies that the di- 
rection of spontaneous magnetization M, parametrized 
by an orientation angle 

= arctan ^ — j , (18) 

begins to change when 6 increases above C , see Fig. 
[U(d). For < # c , the magnetization has only one com- 
ponent M x 7^ with (j> = 0, pointing either parallel 
or anti-parallel to a x which is half-way between cr a {6) 
and cr b (6), see Fig. [TJ Below 6 C the ferromagnetic 
ground state is doubly degenerate and the magnetiza- 
tion is ±M = ±|M X |. When 6 increases above 6 C the 
magnetization begins to rotate in the {M X ,M Z } plane 
by the non-zero angle ±0 Eq. ([T8]) with respect to the 
±|7l/P| initial magnetization below C , and each of these 



two states splits off into two ferromagnetic states rotated 
by ± 1 0| with respect to the a^-axis. As a result, one finds 
four degenerate states above C , and each of them is tilted 
with respect to icr^, either toward ±<r a (#), or toward 
±cr 6 (#), depending on the sign of the rotation angle <\>. In 
the OCM limit 9 = 90° is approached, the magnetization 
angle approaches <\> = 7r/4. In this limit there are four 
degenerate Ising-type ferromagnetic states, with magne- 
tization either along ±a a (90°) (and (a 6 (90°)) = 0), or 
±a b (90°) (and (a a (90 )) = 0). 

Qualitatively the same results were obtained from the 
embedded L x L clusters and they are also shown in Fig. 
[4] for comparison. While 2x2 cluster is too small and 
the quantum fluctuations are severely underestimated, 
the two larger 3x3 and 4x4 clusters are qualitatively 
similar and estimate the QPT point from above, see Fig. 
|4j Rather slow convergence of these results toward the 
MERA result for AE and for \<p\ demonstrates the im- 
portance of longer-range correlations for the correct de- 
scription of the QPT at — C . 

Altogether, these results show that the degenerate 
ground state of the generalized OCM consists of a mani- 
fold of states with broken symmetry. This confirms that 
the OCM is in the Ising universality class, 6,7 with no 
quantum coupling between different broken symmetry 
Ising-type states. However, we found the large value 
of 6 C « 84.8° rather surprising and we investigated it 
further using spin-wave theory. These calculations are 
presented in the next Section. 

Another surprise is the absence of any algebraically de- 
caying spin-spin correlations in the MERA ground state 
at C . They could arise from the subleading eigenvalues 
A2, A3, . . . which we found to be non-zero. However, their 
corresponding coefficients c a p with a > 1 or f3 > 1 in Eq. 
([TT]) are small (at most ~ 10 -4 ) and they decay with in- 
creasing dimension x an d especially the number of non- 
universal layers N. As a result, the only nonvanishing 
term in Eq. ([TT]) is the leading one for a = /3 = 1, de- 
scribing the non-decaying long range order. Notice that 
this observation does not exclude non-trivial short range 
correlations up to a distance 3^ described by the N non- 
universal layers. We believe that when N is too small, 
then the missing short range correlations find a way to 
show up in the small but non-zero universal coefficients 
c a p, but these coefficients decay quickly with increasing 
TV as the short range correlations become accurately de- 
scribed by the increasing number of non-universal layers. 



V. SPIN WAVE EXPANSION 

Since the spin wave expansion in powers of 1/5 be- 
comes exact when the spin S —> 00, we introduced a 
large-^S extension of the generalized OCM Hamiltonian 
Eq. (pQ) with rescaled spin operators: <j x — » S x / S and 
a z — ^ S z /S. We consider first the classical energy per 
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FIG. 5: (Color online) Mechanism of the QPT in the gener- 
alized OCM Eq. (1) for S = 1/2 and = 87° > C . The 
minimum of the classical energy i£o(87°,</>) Eq. ([T9]) (dashed 
line) at cj) — is shallow and thus unstable against weak 
quantum fluctuations which induce two symmetric minima at 
a finite value of ±0 m in obtained from Eq (87° , cj)) derived from 
Eq. ([22])- For better comparison, Eo and Eq are shifted to 
have a minimum value of 0. 

site: 

^(fl,0) = (W(e)) = -i[l + cosflcos(20)] , (19) 

obtained using the mean-field (MF) for the ordered state 
of classical spins 5, with the magnetization direction 
given by Eq. (fl8j) . The classical energy has a minimum at 
<j) = for the entire range of G [0°, 90°). However, when 
the angle approaches 90°, the minimum becomes more 
and more shallow, and finally disappears completely at 
6 = 90°. Thus, the classical ground state becomes very 
sensitive to quantum fluctuations in the vicinity of the 
maximally frustrated interactions in the OCM. 

This behavior of the classical ground state energy ex- 
plains why small energy contributions due to quantum 
fluctuations may play so crucial role in the generalized 
OCM only in the regime of close to 90°, where they 
trigger a QPT by splitting the shallow symmetric classi- 
cal energy minimum at <j) = into two symmetry-broken 
minima at finite values ±(/> m in — we show an example of 
this behavior in Fig. [5] for a particular value of > C . 
Since the quantum fluctuations induce here symmetry 
breaking instead of making the ground state more sym- 
metric, this mechanism goes beyond the Landau func- 
tional paradigm. 

We analyzed the effects of quantum fluctuations 
and the arising symmetry breaking using the Holstein- 
Primakoff representation of spin {5g} operators via {bij} 
bosons: 

cos d> S£ + sin 4> 5?- = S- btb^ , (20) 

&t 

- sin0 Sfj + co S <p Sfj = -fy/2S - b^by + H.c. .(21) 

Operators {b^-jftt-} satisfy standard bosonic commuta- 
tion relations: [bij,bi'j*] = and \bij^b\,^} = 5a'5jj>. In 



this approach, we are looking for a critical value C , above 
which it is energetically favorable to change the direction 
of magnetization M from the symmetric state <fi = to 
a symmetry-broken state with a finite value of <ft ^ 0. 
We expanded the square root in Eq. ([2T]) in powers of 
1/(25) and obtained an expansion of Hamiltonian Eq. 
(PQ) in powers of the operators As we applied 

Wick's theorem to reduce the obtained Hamiltonian to an 
effective quadratic Hamiltonian, the terms proportional 
to the odd powers of 1/(25) do not contribute and are 
skipped below (for more details see the Appendix). When 
truncated at the sixth order term this expansion reads 

H 6 ~H + (2S)~ 1 H 2 + (25)- 2 tf 4 + (2S)~ 3 H 6 . (22) 

Here B.2n is a sum of all terms of the J2n-th order in 
{Kj^lj} operators. In a similar way, and H2 de- 
note expansions truncated at the fourth and second order 
terms, respectively. We have found a posteriori that the 
second order expansion H2 (noninteracting spin waves) 
does not suffice and higher order terms are necessary. 
Consequently, we consider below Hamiltonian Eq. (pQ) 
expanded up to the sixth order. 

For given and 0, we can approximate the ground 
state of the boson Hamiltonian given by Eq. ([22]) by 
a Bogoliubov vacuum obtained as the ground state of 
the quadratic Hamiltonian obtained using the MF 

averaging of four- and six-boson terms. Details of this 
calculation can be found in the Appendix B 

First we performed separate calculations for H 2l H4 
and Hq for several values of spin 5 > 1 when the 1 / (25)- 
expansion given in Eq. ([22]) is convergent. The quadratic 
H2 fails for large where its spectrum becomes gapless 
and the magnetization M Eq. ([T7]) diverges. In contrast, 
and Hq give only small reduction in M in the entire 
range of see Fig. [61(a) and(6^c). Interestingly, the Bo- 
goliubov spectrum remains gapful at C in both the fourth 
and sixth order expansions and, just like in the MERA, 
there are no algebraically decaying spin-spin correlations. 
The critical angle C at which the symmetry-breaking 
QPT occurs increases toward 90° with increasing 5 when 
the quantum fluctuations become less significant. There- 
fore, the magnetization M increases with increasing 5 
and it tends to 1 in the classical limit 5 — >• 00. 

Encouraged by these results, we also performed sim- 
ilar calculations for the generalized OCM Eq. (pQ) with 
5 = 1/2, where the convergence of the l/(25)-expansion 
becomes problematic. Unlike for 5 > 1, we find that the 
fourth order expansion is insufficient as it predicts the 
first order QPT [Fig. E^d)] and does not agree qualita- 
tively with the prediction of the MERA, see Sec. HV1 
Only in the sixth order one finds a qualitative agreement 
between the present boson expansion and the MERA, 
both giving the second order QPT at C . A cusp in M(6) 
seen in Fig. [6](c) shows that even the sixth order ex- 
pansion is not quite converged for 5 = 1/2. Again, the 
Bogoliubov spectrum remains gapful at C in the sixth 
order expansion, with a finite gap equal 1.52, and one 
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FIG. 6: (Color online) Symmetry breaking in the ground state 
as obtained from the boson expansion Eq. (|22|) . Panels (a) 
and (b) show results for S = 1, (c) and (d) — for S = 1/2; (a) 
and (c) depict magnetization M Eq. (|17p , (b) and (d) — the 
value of the magnetization angle <j> Eq. (|18[) that minimizes 
energy. Calculations for He predict the following values of Q c \ 
85.89°, 86.9°, 88.2°, and 89.2° for S = 1/2, S = 1, S = 2, and 
S = 5, respectively (the last two not shown), and C — »• 90° 
for S — > oo. 



finds no algebraically decaying spin-spin correlations. 



Fig. [5l The minimum becomes more and more shallow 
as the compass point = 90° is approached. However, 
the quantum fluctuations are weak due to the gapful or- 
bital wave excitations, and only very close to the above 
OCM point become strong enough to split the shallow 
minimum into two distinct minima in the vicinity of the 
OCM point. In this way the symmetry between the axes 
a and b is spontaneously broken. For this reason the or- 
bital e g model with ferro-orbital interactions, considered 
in Ref. [3l| and corresponding to a "moderate" value of 
— 60° [see Fig. DJb)], orders in a symmetric (uniform) 
phase induced by the stronger (here ex crfjCrfrjr) interac- 
tion component. 

Interestingly, since — unlike in the Landau paradigm 
— the symmetry in the present model Eq. (pQ) is broken 
rather than restored by quantum fluctuations, we do not 
find any algebraically decaying spin-spin correlations at 
the critical point found in the generalized orbital compass 
model Eq. (pQ). The spin waves also remain gapful at this 
point. 
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VI. CONCLUSIONS 



Summarizing, we found that a second order quantum 
phase transition in the generalized orbital compass model 
Eq. (pp) occurs at C = 84.8° which is surprisingly close 
to the compass point = 90°, i.e., only when the in- 
teractions are sufficiently strongly frustrated. There is 
spontaneous ferromagnetic magnetization at any angle 
G [0°,90°]. Below C the ferromagnetic ground state 
is doubly degenerate with the spontaneous magnetiza- 
tion, either parallel or antiparallel to the average direc- 
tion afj + o\y None of the directions, neither a nor 6, 
is preferred in this symmetric phase. In contrast, when 
increases above C the symmetry between a and b be- 
comes spontaneously broken and the ferromagnetic mag- 
netization begins to align parallel/antiparallel to either 
afj or The ground state is fourfold degenerate in 
this symmetry-broken phase. The spontaneous magneti- 
zation M is close to 1 and quantum fluctuations remain 
small in the whole range of 6 G (0°, 90°]. 

These results were obtained using the MERA and the 
mechanism of the QPT was explained within the spin- 
wave theory. For classical spins the minimum of energy 
is at one of the two symmetric states with the magne- 
tization either parallel or antiparallel to af- + of-, see 



Appendix: Details of the spin wave calculation 

In this section we present the details of the spin wave 
calculation. We consider the general case of an L x L 
square lattice, with L being odd for convenience. The 
results presented in Section [V] are obtained after taking 
the thermodynamic limit L — >> oo. 

For given and 0, the ground state of the boson Hamil- 
tonian Eq. ([22]) is approximated by a Bogoliubov vac- 
uum obtained as the ground state of a mean-field (MF) 
quadratic Hamiltonian (to be derived later on). 

Terms #2, #4 and Hq in Eq. ([22]) are given by: 

H 2 = 4[l + cos0cos(20)]^&Jb r 

r 

- Sin 2 {^) - ^ (bj&r+ex + Mr+e x + H.C.) 
- Sin 2 U + ^ J 0^r+e y + Mr+e y + H.C.) , (A.l) 
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Hi = - 4 cos 2 f <p - - ) J2 bibt +e Jrb r + ex 

\ / r 

- 4 cos 2 ^ + ^ ^^ +ey 6 r 6 r+ey 

+ \ Sin 2 " f ) £ {^r±e x + &LeJ + H.C.} 
+ 1 Sin 2 U+l)Y, {^r±e y + ^±e y ) + H.C.} , 

(A.2) 

H 6 = 1 Sin 2 (<A-^)E{ ( 5 ^r) 2 &r(6r±e x + &Lj 

- 26t6i +ex 6 2 (6 r+ex + ?4 +e J& r+ e x + H.c.} 
+ i sin 2 (^ + ^E{ ( 6 r (6 r±ey + bt ±ey ) 

- 2&J&t+e y i£(& r+ e y + &i +ey ) & r+e y + H.C.} , (A.3) 

where r = (i, j), e x = (1,0) and e y = (0, 1). The ± signs 
mean here that both terms, with + and — sign separately, 
must be taken into account. 

To derive the quadratic approximation , we re- 
place the boson terms in H4 and Hq with two-boson 
terms and proper averages by means of the MF approx- 
imation and Wick's theorem. This justifies a posteriori 
why the expansion ([22]) is limited only to the terms with 
even number of boson operators. As an example of this 
approximation, consider one of the contributions to H4 in 
Eq. (|A.2|) : 6j^6 r+ex , which is replaced with a quadratic 
term: 

6j6?6 P+ e x - 2(6 r 6 r ) Mr+e x +2(M>r+eJ b% 

+ <6j& P+ e x } h\ + (b 2 r ) btb r+ex - (btb 2 r b r+Gx ) . (A.4) 



The above replacement procedure leads to six MF pa- 
rameters {mj = {mi, m2, . . . , mo] that should satisfy 
self-consistency conditions. These are in fact all possible 
combinations of operators defined on nearest-neighbor 
sites that cannot be derived one from another by com- 
mutation relations and translational invariance of the lat- 
tice, i.e., mi = (fcjbr), rn 2 = {blb r + Gx ), m 3 = (&t& r+ey ), 
m 4 = {bl),m b = (M r+ e x ), and m 6 (Mr+e y ). 

The obtained Hamiltonian is diagonalized by 

the Fourier transformation followed by the Bogoliubov 
transformation. Fourier transformation which is consis- 
tent with periodic boundary conditions &L+i,j = and 
fri,L+i = h y i has the following form: 

b r = ^$> k e ikr , (A.5) 

k 

where k = (k x ,k y ) is the momentum. In the sum, mo- 
mentum components k x and k y take the values (for odd 
L considered here): 

2tt 2tt L - 1 2tt 

= ±1 , ± — • — . (A.6) 

Diagonalization of H^ F is completed by the Bogoliubov 
transformation: 

6k = ^k7k + ^*k7l k ' ( A - 7 ) 

where the modes and are normalized such that 
|^k| 2 — |^k| 2 = 1- The obtained modes are used to 
calculate new values of the MF parameters {mj (i = 
1,2,..., 6). For instance, one of them reads: rri2 = 
(fyffr r +e x ) = l~2 Xlk \ v k\ 2 cos k x . Starting from random 
values, the above steps are iteratively applied until full 
convergence of all {mj is reached, which results in sat- 
isfying the self-consistency conditions. 
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